In [1]:
set.seed(999)
options(scipen = 9)
options(warn = -1)
source("./environment/libraries.R")
knitr::opts_chunk$set(fig.height = 12, fig.width = 18, fig.dpi = 300)
knitr::opts_chunk$set(warning = FALSE)
All packages loaded successfully
In [2]:
name <- "Kenya_E1"
dataset <- read.csv(paste0("./test/", name, "_processed.csv"))
dataset_c <- data.frame(dataset[7:ncol(dataset)], row.names = dataset$Serial)
proj_coord <- data.frame("Easting" = dataset$Easting, "Northing" = dataset$Northing)
xy_coord <- data.frame("Longitude" = dataset$Longitude, "Latitude" = dataset$Latitude)
dataset_c_closed <- cbind(dataset_c, "Res" = 100 - rowSums(dataset_c))
head(dataset_c_closed)
dataset_spdf <- cbind(proj_coord, dataset_c_closed)
rownames(dataset_spdf) <- rownames(dataset_c)
structures_geojson_path <- file.path("./data", paste0(name, "_topo_lines.geojson"))
structures <- st_read(structures_geojson_path, quiet = TRUE)
structures_utm <- st_transform(structures, crs = 32637) # Transform to UTM Zone 35S (WGS84)
output_dirs <- "./test/ck/"
if (!dir.exists(output_dirs)) {
dir.create(output_dirs, recursive = TRUE)
}
| P | Ca | Ti | Mn | Fe | Ni | Zn | Rb | Sr | Y | Zr | Ba | Mg | Al | Si | K | Res | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | |
| 1 | 0.2731 | 3.4795 | 0.3309 | 0.0603 | 2.7248 | 0.0021 | 0.0033 | 0.0035 | 0.0808 | 0.0017 | 0.012300000 | 0.0933 | 1.227527 | 5.454678 | 21.96007 | 1.971416 | 62.32071 |
| 2 | 0.3835 | 3.9587 | 0.2828 | 0.0566 | 2.4054 | 0.0017 | 0.0031 | 0.0037 | 0.0839 | 0.0016 | 0.004100000 | 0.0955 | 1.234100 | 4.319073 | 18.23019 | 2.227139 | 66.70890 |
| 3 | 0.4158 | 4.0720 | 0.2265 | 0.0617 | 2.1634 | 0.0011 | 0.0033 | 0.0034 | 0.0921 | 0.0012 | 0.004547633 | 0.0755 | 1.316831 | 4.052775 | 17.19447 | 2.489834 | 67.82554 |
| 4 | 0.3704 | 3.9589 | 0.2902 | 0.0601 | 2.6374 | 0.0022 | 0.0035 | 0.0030 | 0.0837 | 0.0021 | 0.010200000 | 0.0880 | 1.384247 | 4.668478 | 18.72138 | 2.425924 | 65.29028 |
| 5 | 0.2981 | 3.3979 | 0.2984 | 0.0707 | 2.7003 | 0.0026 | 0.0030 | 0.0035 | 0.0797 | 0.0017 | 0.015200000 | 0.1005 | 1.280048 | 4.962708 | 20.89937 | 2.486514 | 63.39976 |
| 6 | 0.4373 | 4.0080 | 0.1808 | 0.0537 | 1.9270 | 0.0010 | 0.0029 | 0.0030 | 0.0860 | 0.0012 | 0.003527636 | 0.0583 | 1.304530 | 3.233936 | 13.66489 | 2.973143 | 72.06078 |
In [3]:
source("./utils/functions/create_quick_map.R")
options(repr.plot.width = 10, repr.plot.height = 10)
dbscan_result <- dbscan(dataset_spdf[, c("Easting", "Northing")], eps = 10, minPts = 5)
Area <- factor(dbscan_result$cluster)
dataset$Area <- Area
create_quick_map(dataset, structures, group_data = "Area")
dataset_spdf$Area <- Area
dataset_spdf <- dataset_spdf[dataset_spdf$Area == 1, ]
In [4]:
create_spatial_objects <- function(data_spdf, transformation = "none") {
coords <- data_spdf[, c("Easting", "Northing")]
crs_string <- "+proj=utm +zone=37 +north +datum=WGS84"
if (transformation == "ilr") {
spatial_data <- as.data.frame(ilr(data_spdf[, -c(1, 2, ncol(data_spdf))]))
} else if (transformation == "clr") {
spatial_data <- as.data.frame(clr(data_spdf[, -c(1, 2, ncol(data_spdf))]))
} else {
spatial_data <- data_spdf[, -c(1, 2, ncol(data_spdf))]
}
return(SpatialPointsDataFrame(coords = coords, data = spatial_data,
proj4string = CRS(crs_string)))
}
spdf_comp <- create_spatial_objects(dataset_spdf, "none")
spdf_ilr <- create_spatial_objects(dataset_spdf, "ilr")
spdf_clr <- create_spatial_objects(dataset_spdf, "clr")
source("./utils/functions/lag_distance_from_spdf.R")
source("./utils/functions/site_diagonal_from_spdf.R")
lag_dist <- lag_distance_from_spdf(spdf_ilr)
site_diag <- site_diagonal_from_spdf(spdf_ilr)
In [5]:
source("./utils/functions/geostatistical_workflow.R")
In [6]:
results_ilr_omni <- geostatistical_workflow(spdf_ilr, "ILR_Omnidirectional", lag_dist, site_diag, directional = FALSE)
Selected models:
psill range kappa n
Nug 0.00105033 0.000000 0 136
Gau 0.02873073 3.639113 0 78
Exp 0.01735112 3.593177 0 58
Linear Model of Coregionalization found. Good.
[using ordinary cokriging]
Saved 17 maps to: ./test/ck/ilr-omnidirectional/maps/
V1 V2 V3 V4 V5
ME -0.00002187111 0.02407150016 0.00281125028 -0.00226650633 -0.00243363142
MSE 0.04929753 0.09519844 0.02218845 0.02013800 0.03609853
V6 V7 V8 V9 V10
ME 0.00049795368 0.00378105592 -0.00102092507 0.00646904283 0.00208785007
MSE 0.01606958 0.01892075 0.01062964 0.02852003 0.20423875
V11 V12 V13 V14 V15
ME 0.00311281970 -0.00211847017 0.00763336337 0.00752284718 -0.00555857358
MSE 0.02407379 0.01572140 0.02125674 0.02344569 0.08470890
V16
ME -0.00015847540
MSE 0.01371514
In [7]:
results_ilr_dir <- geostatistical_workflow(spdf_ilr, "ILR_Directional", lag_dist, site_diag, directional = FALSE)
Selected models:
psill range kappa n
Nug 0.00105033 0.000000 0 136
Gau 0.02873073 3.639113 0 78
Exp 0.01735112 3.593177 0 58
Linear Model of Coregionalization found. Good.
[using ordinary cokriging]
Saved 17 maps to: ./test/ck/ilr-directional/maps/
V1 V2 V3 V4 V5
ME -0.0143490895 -0.0104323151 -0.0081959822 -0.0083762146 -0.0182540391
MSE 0.048068220 0.087240677 0.020802082 0.019074672 0.032819468
V6 V7 V8 V9 V10
ME 0.0020152234 -0.0007858862 0.0008997367 -0.0148447382 -0.0304175161
MSE 0.014418483 0.018582763 0.009193277 0.026394854 0.213054365
V11 V12 V13 V14 V15
ME -0.0086429833 0.0055115959 -0.0094178968 -0.0103814351 0.0244454560
MSE 0.021329791 0.014623150 0.023473660 0.026800473 0.082047630
V16
ME 0.0049010919
MSE 0.011445260
In [8]:
par(bg = "white")
options(repr.plot.width = 15, repr.plot.height = 12)
# Compute MAF analysis
g_temp <- create_gstat_from_spdf(spdf_ilr, method = "ordinary")
v_omni <- variogram(g_temp, width = lag_dist / 2, cutoff = site_diag / 3, cross = TRUE)
maf_result <- Maf(spdf_ilr@data, vg = v_omni)
source("./utils/functions/quick_pca_screeplot.R")
quick_pca_screeplot(maf_result)
maf_dataset <- maf_result$scores[, 1:9]
spdf_maf <- SpatialPointsDataFrame(coords = dataset_spdf[, c("Easting", "Northing")],
data = as.data.frame(maf_dataset),
proj4string = CRS("+proj=utm +zone=37 +north +datum=WGS84"))
results_maf_omni <- geostatistical_workflow(spdf_maf, "MAF_Omnidirectional", lag_dist, site_diag,
directional = FALSE, maf_result = maf_result,
resolution = 2)
Selected models:
psill range kappa n
Nug 0.0000000 0.0000000 0 45
Exp 1.0927552 0.9507579 0 34
Gau 0.1560802 9.9853154 0 11
Linear Model of Coregionalization found. Good.
[using ordinary cokriging]
Saved 17 maps to: ./test/ck/maf-omnidirectional/maps/
maf1 maf2 maf3 maf4 maf5
ME 0.003786155 -0.023257962 0.014695922 -0.040568633 0.000701461
MSE 1.0337736 0.9979298 1.0528312 1.3033154 1.3210178
maf6 maf7 maf8 maf9
ME -0.005600492 0.013122820 -0.038253927 -0.031507740
MSE 1.2730757 1.1258029 0.8005447 1.2387432
In [9]:
g_temp2 <- create_gstat_from_spdf(spdf_ilr, method = "ordinary")
v_dir <- variogram(g_temp2, width = lag_dist / 2, cutoff = site_diag / 3,
alpha = c(0, 45, 90, 135), tol.hor = 22.5, cross = FALSE)
maf_result2 <- Maf(spdf_ilr@data, vg = v_dir)
quick_pca_screeplot(maf_result2)
maf_dataset2 <- maf_result2$scores[, 1:10]
spdf_maf2 <- SpatialPointsDataFrame(coords = dataset_spdf[, c("Easting", "Northing")],
data = as.data.frame(maf_dataset2),
proj4string = CRS("+proj=utm +zone=37 +north +datum=WGS84"))
results_maf_dir <- geostatistical_workflow(spdf_maf2, "MAF_Directional", lag_dist, site_diag,
directional = TRUE, maf_result = maf_result2,
resolution = 2)
Directional approach not fully implemented yet - using omnidirectional as fallback
Selected models:
psill range kappa n
Nug 0.9821561 0.000000 0 55
Exp 0.1721345 2.481358 0 35
Gau 0.2113754 13.329098 0 12
Sph 0.4112454 2.030358 0 8
Linear Model of Coregionalization found. Good.
[using ordinary cokriging]
Saved 17 maps to: ./test/ck/maf-directional/maps/
maf1 maf2 maf3 maf4 maf5
ME -0.006610265 0.012043116 -0.034670502 -0.006719259 -0.018199105
MSE 1.0824425 1.1081471 1.0520101 0.9895489 1.1084641
maf6 maf7 maf8 maf9 maf10
ME 0.038734058 -0.048576979 -0.016445294 0.082279189 -0.047689524
MSE 1.0972327 1.1123020 0.7993823 1.1104866 0.9508502
In [10]:
source("./utils/functions/create_ck_maplist_from_list.R")
ck_files <- list.files("./test/ck/", full.names = TRUE, pattern = paste0(name, ".*\\.rds$"), recursive = TRUE)
ck_list <- lapply(ck_files, readRDS)
if (length(ck_list) > 0) {
cat("Available cokriging results:", length(ck_list), "\n")
cat("Files:", basename(ck_files), "\n")
if (length(ck_list) > 1) {
cat("Multiple methods completed. Individual results saved in method subdirectories\n")
}
} else {
cat("No cokriging results found\n")
}
cat("Geostatistical workflow completed successfully!\n")
Available cokriging results: 4 Files: Kenya_E1_CK1.rds Kenya_E1_CK1.rds Kenya_E1_CK1.rds Kenya_E1_CK1.rds Multiple methods completed. Individual results saved in method subdirectories Geostatistical workflow completed successfully!